On the Munn-Silbey approach to polaron transport with off-diagonal coupling and 
temperature-dependent canonical transformations 
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Improved results using a method similar to the Munn-Silbey approach have been obtained on 
the temperature dependence of transport properties of an extended Holstein model incorporating 
simultaneous diagonal and off-diagonal exciton-phonon coupling. The Hamiltonian is partially di- 
agonalized by a canonical transformation, and optimal transformation coefficients are determined 
in a self-consistent manner. Calculated transport properties exhibit substantial corrections on those 
obtained previously by Munn and Silbey for a wide range of temperatures thanks to a numeri- 
cally exact evaluation and an added momentum-dependence of the transformation matrix. Results 
on the diffusion coefficient in the moderate and weak coupling regime show distinct band-like and 
hopping-like transport features as a function of temperature. 

PACS numbers: 



I. INTRODUCTION 

Investigations on organic molecular crystals pioneered 
by Pope and cowokers 1 - more than half a century ago 
made an immense impact on various fields of relevance in 
physics, chemistry, and materials science. The interplay 
between geometric and electronic structures has opened 
up novel materials with unlimited possibilities. The dis- 
covery of conducting properties of doped polyacetylene^ 
in 1977 led to a Nobel prize in chemistry establishing 
a whole new field of organic electronics. After decades 
of research, many families of the organic semiconduct- 
ing materials have been investigated. Among them, 
oligoacenes*£~~ and oligothiophenes^ have been examined 
intensively due to the possibility of making single crys- 
tals with few defects via repeated sublimation processes. 
Single crystals of these materials can be used to develop 
transistors that allow for fundamental research on intrin- 
sic properties as well as performance of organic electronic 
devices. The Holstein molecular crystal model was intro- 
duced in the 1950sZ£ to account for the effect of electron- 
phonon interactions on transport. Since then, a lot of 
work has been carried out to deal with the complexi- 
ties of the transport phenomenon in organic molecular 
crystals. As reviewed in ReL 1 ^, charge transport in or- 
ganic semiconductors is governed by electronic coupling 
and electron-phonon interactions, and the main trans- 
port mechanism can be described by polaron and dis- 
order models. Many theories on polaron transport are 
based upon a phenomenological model^ or utilize the po- 
laron effective mass approaches. However, validity of 
these models is restricted to specific parameter ranges. 
More general microscopic models in understanding po- 
laron transport was develop ecUi~— from a density ma- 
trix approach, which is capable of describing electronic 
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coupling, diagonal and off-diagonal electron-phonon in- 
teraction of arbitrary strength over a wide range of tem- 
peratures. The effect of off-diagonal coupling poses much 
more difficulties to tackle than its diagonal counterpart, 
and as a result, off-diagonal coupling has not been exten- 
sively treated before 14 . More recently, microscopic mod- 
els from generalized master equation approach 15 , and dy- 
namical mean-field theory 1 - have been developed. 

Polaron transport theory developed by Holstein in his 
seminal worki£ was based upon perturbation expansions 
with its validity sometimes limited to cases of small trans- 
fer integrals^. Despite its limitation, the theory has been 
widely used for qualitative interpretations of experimen- 
tal data, including the temperature dependence of band- 
narrowing effect as well as the crossover from band-like 
behavior to hopping transport with increasing tempera- 
ture. A microscopic transport theory that accounts for 
simultaneous diagonal and off-diagonal electron-phonon 
coupling was developed by Munn and Silbej^ 1 - - — based 
on a canonical transformation of the Hamiltonian with an 
element of optimization. Off-diagonal coupling is found 
to increase the polaron binding energy and introduce new 
minima and broadening to the polaron band. In con- 
trast, the effect of diagonal coupling is restricted to band 
narrowing. The Yarkony-Silbey variational approac h 20 i 21 
provides an attractive direction to solve the transport 
problems in organic molecular crystals. This approach 
has been borrowed to treat two phonon bands in a three- 
dimensional lattice by Parris and Kenkre 2 ^. However, 
as the Yarkony-Silbey approach lacks the flexibility in 
variational parameters, a more sophisticated approach is 
needed to capture the correct temperature dependence of 
diffusion coefficient and mobility. To obtain the optimal 
basis for a polaron system, finite-temperature variational 
method by Cheng and Silbej*^ combining Merrificld's 
transformation with Bogoliubov's theorem has been re- 
cently devised bringing substantial improvements to re- 
sults obtained previouslj* 1 ^. The theory is able to capture 
the universal band-like-to-hopping transition as temper- 
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ature increases, and more importantly, a temperature- 
independent mobility at extremely low temperature has 
been found in agreement with a nonperturbative ap- 
proach recently developed by Ortmann et al 22 ^ 23 . How- 
ever, Merrifield's transformation itself is applicable only 
in the small polaron regime, and the off-diagonal exciton- 
phonon coupling was not accounted for in Ref. 18 , making 
the approach insufficient to describe the phonon assisted 
transport /tunnelingi^ in molecular crystals. Thus, fur- 
ther improvements are necessary for cases with simul- 
taneous diagonal and off-diagonal coupling. Following 
up on the Munn-Silbey transformation method, Zhao 
et al^r— devised a self-consistent routine to determine 
the transformation coefficients, demonstrating substan- 
tial corrections of the polaron band structure thanks to 
a built-in momentum dependence of the transformation 
coefficients. This serves as the starting point towards the 
evaluation of transport properties in this paper. 

An attempt to extend the Holstein model to higher 
dimensions with a microscopic model has been made by 
Kenkre et al.,— utilizing the generalized master equa- 
tion approach. This model was able to adequately ex- 
plain the temperature dependence and anisotropy of mea- 
sured mobilities in naphthalene. Using master equations, 
Wang et al.^ developed a nonperturbative method 
to handle the electron-phonon coupling fully quantum- 
mechanically and quantify charge-carrier transport prop- 
erties in organic molecular crystals. A more sophis- 
ticated microscopic model based on a Hamiltonian of 
the Holstein-Peierls type has been presented by Bob- 
bert and colleagues^Sr^i. Very recently, a theory based 
on non-perturbative evaluation of Kubo formula for 
the carrier mobilit y 22 ' 32 ! 33 has been proposed, showing 
several improvements including the elimination of low- 
temperature singularity that often appears in theories 
based on narrow-band approximations. Based on a non- 
perturbative evaluation of the Kubo formula, calculations 
have been carried out for durene crystals 34 and guanine 
based materials 3 -^, and transport channels of charge car- 
riers are revealed by making use of ab-initio tools. Cou- 
pling to the intermolecular acoustic modes^ was found to 
play a significant role in charge transport, thus the model 
presented by Bobbert et al. may be improved further. 
Experimental data of naphthalene can be well reproduced 
by this model with microscopic parameters obtained from 
ab-initio calculations. However, this model neglects the 
intra-molecular modes as well as the coupling of excitons 
to the intermolecular acoustic phonons, and only the cou- 
pling to the optical modes is accounted for. Moreover, 
the models developed by Bobbert et al., and by Munn 
and Silbey, are based on nonlocal canonical transforma- 
tions with additional approximations. Thus the range of 
validity of these models still requires further exploration, 
despite some qualitative agreements with experiments. 
Hence, in the light of recent theories by Hultell et a/.— 
and Troisi et alM&, it is clear that both local and nonlo- 
cal exciton-phonon interactions should be taken into ac- 
count. A comprehensive theory of the charge transport in 



organic crystals requires treatments of the Hamiltonian 
having the lattice dealt with quantum mechanically in 
order to make itself valid for any temperature of interest, 
and the Munn-Silbey approach 1 ^— combined with the 
self-consistent routine devised by Zhao et al*^ provides 
a good opportunity to explore the simultaneous effect of 
diagonal and off-diagonal exciton-phonon coupling on the 
transport properties of polaron. 

The rest of the paper is organized as follows. The the- 
oretical models arc firstly introduced in Section II, where 
the essential improvements over the original Munn-Silbey 
model have been proposed. In Section III, the numerical 
results obtained with our approach have been discussed 
in detail and comparisons with previous results are made 
to illustrate the success of our approach. Finally, the 
conclusions are drawn in section IV. 



II. METHODOLOGY 

The generalized Holstein Hamiltonian and its trans- 
formation have been described in detail previously by 
Munn and Silbeyii~— . Here we will first discuss briefly 
the canonical transformation of the Hamiltonian which 
yields a weak residual excitation-phonon coupling. Eval- 
uation of correlation functions, essential to the final cal- 
culation of diffusion coefficient and mobility as a function 
of temperature, follows next. Our major focus in this pa- 
per will be the explicit form of the correlation functions, 
evaluation of critical exponential matrices and especially 
the transformation coefficients without resorting to the 
approximations in Ref* 1 ^. In our study, the transforma- 
tion coefficients are determined through solving a set of 
self-consistency equations, rather than an analytical form 
approximated in Ref* 1 ^, therefore providing substantial 
corrections in the transport calculations. In Sec. II A 
and B, expressions for the scattering and the hopping 
rates as well as the diffusion coefficients are introduced. 
Correlation functions are presented in Sec. II C. 

A. Hamiltonian and its Transformation 

The generalized Holstein Hamiltonian incorporating si- 
multaneous diagonal and off-diagonal coupling— has the 
form: 

H ^ ^ ^a^a n -f- ^ ^ Jmn^n^rn 

n n : m 

+ 5> (? (&t& g + l/2) 

Q 

+N- 1 / 2 Y, u q fl n ala m (b q + bl q ), (1) 

nmq 

here at(a„) is the creation (annihilation) operator of 
an excitation (i.e., exciton or charge carrier) with en- 
ergy e at site n, and b q (b q ) creates (destroys) a phonon 
with frequency uj q and wave vector q. J nm stands for 
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the electronic transfer integral between site n and m. 
In this study, the discrepancies between molecules have 
been neglected, thus the on-site energies are set to e 
as shown in EqfTJ Finally, the excitation-phonon cou- 
pling is described by the last term of Eq|T] where the 
coupling strength /* m must satisfy the relationship : 
(f«_ m )* = e -iq-(Rn-Rm) f-Q^ to ma k e H Hermitian 
and translationally invariant. 

A unitary transformation can be applied to the Hamil- 
tonian, as shown below: 



with 



with 



H ->■ H = U^HU. 



JJ = e N- 1/2 £ nm , K m (bt q -h)aU 



(2) 



(3) 



The transformation requires f/' = —U in order to 
keep the translational symmetry of A q lm , similar to f% m . 
Details of the transformation have been explicitly intro- 
duced in Ref^-i^. 

Moreover, the Hamiltonian in momentum representa- 
tion is more convenient than the site representation due 
to the translational symmetry. Eq. [T] can be written in 
the momentum representation as: 

H = ^2e k ala k + ^2u} q {blb q + 1/2) 

k q 

+N- i / 2 j2^f-A + M b ^ + b -^ ( 4 ) 

kq 



with e k = e + Jfe , and 

° k — / c ° nm i 



(5) 



where J n . m — J(S n ,m+t + ^n,m-i)j an d w g is the phonon 
frequencies. 

The excitation-phonon coupling coefficients can also be 
written as: 



fq _ -ife-(R„-R m ) f q 

Jk / j ° Jn—mi 
m 

with following coupling geometry: 

fl=g- «0[sin(fc) - sin(fc - q)], 



(6) 



(7) 



which indicates that the antisymmetric type of coupling 
is used throughout the paper. 

Now the transformation takes the following form in the 
momentum representation: 



U = e 



(8) 



and 



A l = ]>>- l ( k - q > R "e* k - R ^ m - {A- k l q )*. (9) 



The transformation can then also be written as: 

[/ = e E fe A,'4 s M'«fc' ) (10) 



S kk ,=N-V 2 A k _- k k '(bt _ k -b k . 



(11) 



where Sfcfc' is an operator creating a net phonon with 
momentum k! — k. This operator is essential towards the 
evaluation of a series of important matrices, where 



a k 



^ ^fcfc' a k ' 



with 



6» fcfc / = [exp(-S)]jwfe/ 
4, = [exp(S)W 



(12) 



(13) 



where S stands for matrix S kk > ■ The transformed Hamil- 
tonian can be partitioned into Hq terms (zeroth order 
part) and the perturbative V terms. Thermal average of 
the 6kk' operator has the property: 

(Olk+q^M+q') = V. (I 4 ) 

which ensures the correct thermal equilibrium behavior 
of Hq by keeping it diagonal in excitation wave vector. 
Ho here has following expression: 

H = 5>,(&t&, + 1/2) + J2 

a[a k {(e - J k ) 

q k 
q 

X Y,^f-k' A - q k(0l + q^k,k-q)+N-^ 
qk' 

X J2 U «f-V A -U€+q, k 0k'Abq + tf-q))}> ( 15 ) 
qk' 

here J k represents the renormalized transfer integral 
which can be expressed by: 



Jk = ^2,Jk'{0 yk 0k'k)- 



(16) 



The perturbation parts of the transformed Hamilto- 
nian can be written as: 

V = UkTkk'-kk" - 2N ^ 1 ^2 UJ qf-k( A -k") 



x T k 



k-\-q,k' :k,k' 



a k " 



x [T k +qM'-kk"(b q + bLg) - (T k , 

-\-q,k' -kk" 

+ N-y"Y.^ A -k-Y,f-M' +q ^,k-q)] 

qk k' 



x (b q + bl q )al +q a k , 



(17) 
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where 

Tkk' ;uu> = @ kk i 9uu' — (Qkk>® uu ')- (1®) 

However, the last term of V has the potential to grow 
with increasing temperature, motivating the choice of A k : 

A q _ k = J2f- k M' +g , k ^,k- q ), (19) 
fc' 

to curb this possible uncontrolled growth of 9 k k' ■ 

By applying the thermal average routine shown in 
Refii^, the transformed Hamiltonian can be written in 
a simpler form, which facilitates the numerical calcula- 
tions. The zeroth order Hamiltonian Ho now can be writ- 
ten as: 

Ho = ^hakak +^2^ q (b\b q + -), (20) 

k q 

and 

~e k = e + Jk-N- 1 J2\A q k \ 2 u Jq . (21) 
q 

The renormalized transfer integral in Eq. [21] is given by: 

Jk = ]T(# fe ) 2 (exp£V,J fe ,, (22) 

fc' 

where the triadic E kk , is a quantity introduced to yield 
a simpler form of the transformed Hamiltonian, with fol- 
lowing relationship: 

E\ k , = N-\2n k -k' + l)(4lf TA k k r k \ (23) 

here n q is the Bose-Einstein distribution, with 2n q + 1 = 
coth(±/37iu;). 

With the form of A q k shown in Eq. [T51 the perturbation 
part V of H can be written as: 

V = ^ {JkTkk'-kk" — 2iV 1 

x "if-k A -l» T k+<i,k';k,k»- q + N~ 1/2 

x Y,^-k T k+ q ,k':kk''(b q + bl q )alak''}. (24) 

In this study, the transformation coefficients A k (A- 
matrix) will be obtained variationally prior to the cal- 
culations of the transport properties. The detailed nu- 
merical variational procedures can be found in Rcf. 24 . 
In performing the variational calculations, the transfor- 
mation coefficients must satisfy following self-consistency 
equations together with Eq. [23] 

(e k )=exp[-^E° kk ,), (25) 
fc' 



A\ = (e k - q )(6 k ) £ f q k , {eME q )}kk>- (26) 

fc' 

The self-consistency equations above follow the Munn- 
Silbey secular elimination scheme. In Ref»A2,, this set of 
equations was solved upon further approximates on A k , 
whereas in this paper, it is solved numerically with an 
accuracy unavailable previously. 

In the Munn-Silbey approach^, the transformation co- 
efficients A q k was at first approximated by the scaling pa- 
rameters £ and r\ as: 

A q k =g^-iHMk)-Mk-q)]- (27) 

Numerical evaluation of the Eqs. l2"3"]|2"6"l is facilitated by 
rewriting the Eq. [27] with real matrices £| and r) k in the 
following form: 

A\ = 9% - iW k [sin(fc) - sin(fc - q)] . (28) 

In the equation above, the sine function leads to zero 
imaginary part along the line with q — 0, 2k ± n. In the 
numerical calculation, the lines with q — and q = 2k±ir 
are treated as removable singularities. Moreover, the val- 
ues of r) k along these lines are chosen to be analytically 
connect with neighboring values^. The variational ma- 
trices £j? and vi k chosen can preserve the symmetry prop- 
erties of A k in a way such that £j? = £j~l q and X]\ = T) k q ■ 

The self-consistency approach shown in this paper is 
largely limited by the secular elimination scheme that 
is independent on the rigid-lattice tunneling matrix el- 
ements "J". However, since "J" determines the rigid- 
lattice band structure and related quantities such as effec- 
tive mass and density of states, this approach only applies 
for the narrow (rigid) band regime where the traditional 
small polaron theory is still available. The wide (rigid) 
band systems are beyond the scope of the approach dis- 
cussed in this work. Based on the simple introduction 
and discussion of the Hamiltonian and transformation 
schemes, the parameters used for our numerical calcula- 
tions are limited to the regime which holds the traditional 
small polaron scenario. 

B. Diffusion Coefficient 

The diffusion coefficient is closely related to the mean 
square displacement of an excitation in a given period of 
time. For the transformation scheme shown in this paper, 
the mean square displacement can be calculated from 
the excitation density matrix. Approximations have been 
made to reduce the complexity of the exact formulism of 
the mean diffusion coefficient by neglecting small terms 
of the perturbation V as shown in the section above. The 
diffusion coefficient can then be expressed a o n i 13 : 

# = a 2 ((^/rfcfc+7fcfc», (29) 

where a is the nearest neighbor distance, k is the wave 
vector, and v k — VkE k - The double bracket denotes the 
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thermal average of the polaron states E k , where T kk and 
7^^ are the scattering and hopping rates expressed as 



(30) 



7fcfc 



-v 2 V 



Re ("W^fcjfe'+^fc'+g — W / "fe,fc+g;fc',fc'+ ? )| ? =0-(31) 

where W/cfc+^fc'.fc'+g, crucial to the calculation of diffu- 
sion coefficient, are given by 



W, 



k,k-\~q;k' ,k'-\-q 



dt{(V k > +q .k+ q Vk,k'{t))exv[~i{E k , +q - E k+q )t] 



+ (V k > +q , k+q (t)V kk >)exp[-i(E k - E k ,)t}}, 



(32) 



here the single angular bracket denotes the average of 
the residual interaction V kk i(t) over phonon states. Us- 
ing the formulism above, the evaluation of D will be done 
numerically here without introducing further approxima- 
tions as in Re h 11 ' 12 . With the diffusion coefficient D, the 
mobility can be easily obtained with the Einstein relation 
fx = e/3.D. 



C. Correlation Functions 

As inferred from the expression of the diffusion coef- 
ficient, the central problem is the evaluation of the cor- 
relation functions (V k >+ qyk + q V k . k > (t)) since both the scat- 
tering and the hopping rates are determined by the W 
quantities. It is indicated from the previous formulism 
that the residual interaction V kyk ' (t) is closely related to 
the transformation coefficients A q k . Hence, the problem 
of evaluating D ultimately hinges on the solution of the 
self-consistency equations, i.e., Eqs. |2"31 [231 and 121)1 In this 
paper, the transformation coefficients that determine the 
structure of polaron states and energy bands are solved 
numerically, which provide new information about the 
transport properties of polarons in molecular crystals. 

Without the analytical approximation form of r\ carried 
out in Ref J^, the important quantity (6 k k , 9 qtq > ) can be 
derived asi2 : 

(0iU<V) = (e- k ')exp{D k - q ) k , k '{e- q >}5 k - k , >q - ql 

= cxp(£ fc -«)_ g ,,_g(fl_ g ,)4-fc', g -g', 

(33) 



and 



M = jfPk-k'(t)(A k k Z k q 'yA k k - k ', (34) 
(D%Mt] = ^Pk-k'{t){Atl k 'TA k _- k ' +q , (35) 
where the expression of (9 k ) is given by Eq. 1251 



The factors Pq (t) are the characteristics of correlation 
functions for linear exciton-phonon coupling 39 , which 
may take the form 

P Q (0) = (2n Q + l)/N 

P Q (t) = [{2n Q + l)cos(wt) + iwa(ut)]/N, (36) 

where uq is also the Bose-Einstein distribution, with 
(2uq + 1) = coth(i/3fio;) and ui is the mean phonon fre- 
quency. 

In order for the correlation to decay to zero at long 
times, Pq(£) should be multiplied by a decay factor A(t), 
defined as the Fourier transform of the phonon density 
of states, which takes a Gaussian form^i as given by: 



Ait) = cxp(-A 2 i 2 /4), 



(37) 



where A is the phonon bandwidth taken much less than 
lo. For the convenience of numerical calculations, the 
mean phonon frequency lo is set to 1. For all calcula- 
tions performed in this paper, A is set to be 0.1, indi- 
cating a rather narrow phonon bandit. The inclusion of 
phonon bandwidth is a natural result of exciton-phonon 
coupling. For the Holstein polaron model, it has been 
inferred 8 -^ that phonon dispersion is directly linked to 
lattice relaxation time, and it is also crucial for the cal- 
culation of probabilities of hopping events accompanied 
simultaneous absorption and emission of phonon. It is 
also pointed out by Holsteir* 8 - that replacing the phonon 
dispersion by a single frequency (Einstein model) yields 
meaningless results for probabilities of hopping events. 
A number of factors can result in a finite bandwidth in 
the phonon density of states, and as pointed out by other 
author o 11 ! 18 , the inclusion of the phonon bandwidth en- 
sures the long-time decay of correlation functions. 

In the evaluation of (V k '-\-q lk -\- q V klk '(t)}, an important 
quantity is the four 9 k correlation function which is given 
by: 



(fyfcl.fca ^3 > fe 4 Ql ,12 (*)0«3,<?4 (*)) 

^ (^L+r-i M ® k 3 .*s-ri ) (@ q2 +r 2 ,q 2 ® 1 



[exp(-E 



k 1 +q i +r 3 -q 2 -q 3 



(*))]*3 



+r 3 -qi-k 4 -ri,-q 2 ~r 2 



X [eME^-^m-q.-r^-q, 

x [exp(E k ^-^(t))]^ qiM+r ^ qi ^ ri 

x [ exp (-E*-«- r *(t))]- q *,-«-T a 

x S, 



ki +91 +ki +(J 4 ,k 2 +q 2 +k 3 +q 3 



(38) 



Now we can derive the expression for 
(Vfe'+g ! fc+ g Vfc ! fc'(i)), which is the core item of W. 
Here we first express Vk,k' as follows: 



Vk,k> - '^ J {j K T Kk , Kk : - -j=y^^ q fi K 



X [ y—{A_^.T K ^-q ik - K! k> — q + A q _ k T Ki ^. q _ k ^ q ^ Ktk i ) 



;{T K+q ^ k . K ^ k i%l) q + ip q T K+qtk;K!k /)]} 
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with 

T k ,k>;q,q> = #1,*;' <W ~ (^fc,fc'^?,?')) ( 39 ) 

and Jfc = Jcos(fc) denotes the polaron band with its 
width determined by J. As described in the previous 
section, the currently applied approach is only applicable 
for narrow band rigid lattice system. Consequently, the 
value of J should not be too large in order to guarantee 
the accuracy of the transport results obtained. 
The correlation function is now given by: 

(v klM v quq2 {t)) 

= C^2{ j k 3 T k3ikl . k3ik2 - -j=^2,u Ql f_i 3 

k 3 V Qi 

+ A_ 1 k *T k3+ Q l! k x +Q 1 ;k 3 ,k 2 ) 

~ -^{Tki+QukvMMi'Qi + ^QiTka+QiM-Mfo)]} 

X [^j^{^— < qi r ^k i +Q 2 .qi;k i ,q2-Q2 
+ AS_ qi T ki+ Q 2 ^ 1+ Q 2 . ki! q 2 ) 

,qi; k 4,q2 )]})■ 

(40) 

The explicit expression of the correlation function has 
been given in the Appendix A. 

As read from the expression of the correlation func- 
tion in Eq. [40] matrices A\ play an essential role towards 
an accurate evaluation of the correlation function. In 
RefJ^, these matrices were obtained only approximately, 
while in this paper, they are obtained numerically with 
unprecedented accuracy. One will find that the accu- 
rate evaluation of the A q k has dramatic influence on the 
overall transport calculations in the following section, es- 
pecially for the cases with off-diagonal excitation-phonon 
coupling. 

III. RESULTS AND DISCUSSION 
A. Evaluation of Transformation Coefficients 

The numerical evaluation of transformation coefficients 
A q k involves the solution of the self-consistency equations 
Eqs. B5l and BB1 Details of the evaluation procedures 
have been discussed in detail in Reb^i, here we only re- 
visit some important issues related to the transport prop- 
erties, especially the role of off-diagonal coupling. 

For moderate temperature with k^T = Huj, the self- 
consistency method applied in this paper tends to fail 




FIG. 1: A-matrix parameters f| (a)-(c) and r\\ (d)-(f) as func- 
tions of increasing diagonal coupling strength g 2 . g 2 =0.1 for 
(a) and (d), 0.3 for (b) and (e), and 1.0 for (c) and (f). The 
other parameters are: J=0.1, cj> 2 =0.3 and T=1.0. 



in strong coupling regimes, this is attributed to the fail- 
ure to minimize errors to desired precision or converg- 
ing to multiple solutions that are sensitive to procedure 
initialization. It is also possible for the convergence to 
fail in strong coupling regimes for solely numerical rea- 
sons. This problem increases with increasing excitation- 
phonon coupling, which lead to large positive and nega- 
tive eigenvalues for the exponential matrices, thus to the 
possible failure of obtaining reliable diffusion coefficient. 
Due to the problems of convergence and possible multiple 
solutions of the method applied, we restricted our studies 
in the weak and intermediate coupling regimes to ensure 
the robustness of the results obtained. 

As demonstrated in Refi^i, A-matrix parameters 
and rff, in Eq. [25] obtained with the self-consistency pro- 
cedure have a significant wave vector dependence, which 
is lacking in the original treatment of Munn and Silbeyi^ 
as indicated in Eq. [57] In addition, the Munn-Silbey 
parametrization of the A matrix, r\ and £, can be regarded 
as the lower bounds of £| and r\\, rather than the average 
values. More importantly, the wave vector dependence of 
the variational parameters lead to significant changes in 
the structure of polaron energy band E k , where the bi- 
modal variation of E k can be always found. Such changes 
in the structure of polaron states and the energy band 
lead to critical modifications of the transport properties, 
which will be explicitly examined in the next section. In 
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FIG. 2: Polaron bandwidths (a) as functions of tempera- 
ture, and band structures (b)-(c) comparisons between this 
work and those from Toyozawa variational method with off- 
diagonal coupling strength (b) </> 2 =0.1 and (c) (f> 2 —0.3 at T=0, 
where g 2 and J are both equal to 0.1 for all calculations. The 
inset of (a) is the band structure comparison between the re- 
sult from Ref.— by Munn and Silbey and that from this work 
at T=l with g 2 = 2 = J = 0.1. 



this section, we present detailed structures of A-matrix 
parameters as well as the corresponding polaron bands 
obtained with the self-consistency routines, in order to 
facilitate understanding of calculated transport proper- 
ties, which are deferred to the next subsection. 

As shown in Fig. [TJ the two sets of A-matrix param- 
eters in Eq. 1281 K and rfl, differ in their wave vector 
dependence. In the vicinity of q = 0, £| is weakly struc- 




FIG. 3: A-matrix parameters £| (a)-(c) and 77J? (d)-(f) as func- 
tions of increasing temperature. T = 0.2 for (a) and (d) , 1.3 
for (b) and (e), and 4.0 for (c) and (f). Other parameters are 
J=0.1, g 2 =0.1 and <£ 2 =0.3. 



tured and relatively flat with respect to k, while the wave 
vector dependence of 77^ is much more pronounced than 
that of £?. Regardless of coupling strength, the wave vec- 
tor dependence of and has a characteristic "manta 
ray" shape, a fact that may be understood through an ex- 
amination of the real-space structure of the A matrix 24 . 
As revealed in Fig. [1] the A-matrix parameters tend to 
have less wave vector dependence for larger g. If diagonal 
coupling strength exceeds that of off-diagonal coupling, 
the parameters lose most of the q dependence. Thus, the 
transformation coefficients in Eq. 1271 as approximated by 
Munn and Silbeyi 2 -, represent the A-matrix in Eq. [28] in 
the limit of strong diagonal coupling. It is also interest- 
ing to look into the change of the polaron band structure 
(especially, band narrowing) as a function of the tem- 
perature, which may provide an intuitive understanding 
of the difference between transport properties obtained 
here and those by Munn and Silbey using Eq. [27l Vari- 
ance of the polaron bandwidth as well as the entire band 
with respect to temperature has been obtained with our 
numerically determined A-matrix, and is illustrated in 
Fig. O It is also indicated in the inset of Fig. |2[a), the 
polaron band obtained with the transformation coeffi- 
cients of Munn and Silbeyi 2 - lacks a bimodal structure, 
and has a much larger bandwidth than that from the cur- 
rent approach. It is indeed counterintuitive to see that 
the greater momentum-space variability we find in our 
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transformation coefficients should translate into a nar- 
rower polaron band. The explanation lies in that the 
binding energy is proportional to the average of \Al\ 2 . 
Compared with results from Ref.— , the average binding 
energy obtained in the current approach is significantly 
larger and the average Debye- Waller factor significantly 
smaller, which eventually leads to the narrower band. 

To lend further support to the polaron band structure 
calculated here, it is helpful to introduce the Toyozawa 
Ansat a 41 ' 42 for zero-temperature band comparisons. As 
shown in Figs. [2] (b) and (c), good agreement between 
the current method and the Toyozawa Ansatz has been 
found. For weak diagonal and off-diagonal coupling as 
is the case in Fig. [5] (b), the Toyozawa method yields 
a larger bandwidth, and with increasing (j), the bimodal 
structure appears and the agreement between the two ap- 
proaches improves. The concurrence provides legitimacy 
to the current approach in finding reliable polaron bands, 
a step that is essential for the transport calculations. 
Moreover, the band narrowing effect is in good quali- 
tative agreement with those obtained with the Holstein- 
Peierls model and ab-initio calculations in Ref<2£. As 
can be readily seen from Fig. [2] (a) , the off-diagonal cou- 
pling strength helps boost the bandwidth for tempera- 
tures lower than u, g 2 — 0.1, and J = 0.1. 

To better understand the bandwidth narrowing effect, 
the parameter matrices and rfj. have been plotted in 
Fig. [3] for three temperatures, T = 0.2,1.3, and 4.0 in 
unit of u>. Other control parameters in Fig. [3] are set 
to be the same as those for Fig. EJc). As temperature 
increases, the wave vector variations grow for both ^ 
and rjl, while in Fig. [TJ the variations are reduced as 
g increases. As explained in Ref.^ 4 -, a larger wave vec- 
tor variation in £? and "qj. is transformed into a weaker 
distortion of the energy band, which explains the temper- 
ature dependence in Fig. [2] (a) . In contrast to an overall 
enhancement of polaron band-narrowing by including off- 
diagonal coupling in Ref. 2 ^, our model demonstrate the 
possibility that at low temperatures off-diagonal coupling 
may increase the bandwidth while at high temperature 
the effect is reversed. Such an effect can be understood 
from Eqs. |2"D|21)1 where the bandwidth is shown to be 
closely related to the renormalized transfer integral Jk 
and the polaron binding energy, which is the third term 
on the right hand side of Eq. [5T] At lower temperatures, 
the binding energy dominates in Eq. [5T] in the presence 
of off-diagonal coupling strength </>, resulting in a band- 
width that increases with <j>. At high temperatures, how- 
ever, the renormalized transfer integral Jk dominates at 
weak off-diagonal coupling, and the bandwidth derived 
from the Jk term shrinks with increasing <fi. Stronger 
off-diagonal coupling leads to stronger correlation which 
helps to form more spatial extended polaron structure 2 ^, 
which naturally leads to higher coherent phonon-assisted 
exicton transfer, thus the polaron bandwidth at low tem- 
perature for stronger <fi case is lager than the smaller 
one. However as temperature increases, the inclusion of 
stronger off-diagonal coupling strength helps to increase 



the thermal fluctuation of exciton (or polaron) through 
the renormalized J term in Eq. 1211 thus in this case 
phonon has more significant role in reducing the band- 
width. 



B. Transport Properties 

The influence of the diagonal coupling strength g and 
off-diagonal coupling strength 4> (especially the latter) 
is substantial in determining momentum space modula- 
tions of transformation coefficients At , which in turn de- 
termine the structure of polaron states and the energy 
bands. For all calculations here, the number of site is 
equal to 6 and the hopping integral J is 0.1 in unit of 
oj unless specified otherwise. Comparisons between the 
results obtained here and those from Ref^ have been 
made to illustrate the importance of the A-matrix mo- 
mentum dependence as well as the straightforward nu- 
merical evaluation of transport properties sans additional 
approximations. 

In the calculations of the transport properties, com- 
putation time increases dramatically with the number of 
sites. Our existing computational capabilities allow only 
calculations for less than 10 sites without overshooting 
memory requirements. To ensure the robustness of the 
results obtained with current method, we computed the 
diffusion coefficient D for ID lattices of 4, 6 and 8 sites. 
The results with various lattice sizes and control param- 
eters are displayed in Fig. 2J indicating that almost iden- 
tical results are obtained for lattices of 6 and 8 sites in 
the absence of off-diagonal coupling. Hence, in such a 
scenario, calculations with 6 sites are expected to give 
sufficiently reliable values of D at affordable computa- 
tional cost. As (j) 2 increases, however, the number of 
sites becomes an issue in determining the value of D due 
to excitation-phonon correlation spans over longer dis- 
tances. The effect of the lattice size is more prominent 
in the low temperature regime as shown in Figs. EJb) 
and (c), where the transport properties are dominated by 
coherent or band-like contributions. At higher temper- 
atures, the hopping contribution eventually dominates 
the transport and the results become less sensitive to the 
system size. This can also be understood from polaron 
band narrowing at high temperatures and a reduced ef- 
fective transfer integral with increasing temperature. As 
the diagonal coupling strength g increases, the difference 
between the results obtained with 4 sites and 6 sites be- 
come less obvious than the weak coupling cases. Even 
better agreement is reached for <j> 2 = 0, and the results 
are found to be almost identical from 4 to 8 sites as shown 
in Fig. |U[d). It is also interesting to spot nonmonotonic 
changes of the diffusion coefficient (i.e., a "hump") as a 
function of the temperature, in the absence of any size 
dependence, in the low temperature regime ioig 2 = 0.5 
and nonzero c/> 2 , as shown in Figs. 0] (e) and (f). Such a 
hump also appears in Fig.@](c) for g 2 = 0.1 and <fi 2 = 0.7 
despite a slightly visible size dependence of the diffusion 
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FIG. 4: Diffusion coefficient D versus scaled temperature 
&bT/oj for a 1 dimensional chain with 4, 6 and 8 sites, where 
the diagonal coupling strength g 2 = 0.1 for (a)-(c), and g 2 = 
0.5 for (d)-(f). The off-diagonal coupling strength <f> 2 is equal 
to 0.0 for (a) and (d), 0.3 for (b) and (e), and 0.7 for (c) and 
(f). Finally, the transfer integral J=0.1 for all cases. 



coefficient. For a fixed value of g 2 , the "hump" becomes 
more pronounced with increasing 2 , underlying the par- 
tial role of off-diagonal coupling. 

As shown in Fig. [SJ increasing the phonon bandwidth 
A leads to an increase in the band-like contribution to the 
diffusion coefficient, and a decrease of the relative impor- 
tance of the hopping contribution. These findings from 
our numerical calculations for a case of diagonal exciton- 
phonon coupling are in good agreement with analytical 
results from RefJ^. Qualitatively similar results can be 
obtained with the inclusion of off-diagonal coupling. 

In the original work of Munn and Silbey, the diffusion 
coefficient as a function of temperature have been stud- 
ied with structureless A-matrix parameters £ and r\. The 
detailed structure for the A-matrix parameters is found 
to have a substantial impact on the transport properties, 
as revealed in Fig. [5] by comparing the diffusion coeffi- 
cients from Ref.— with those obtained here. Moreover, 
we would like to note that in the absence of off-diagonal 
coupling our results are in good agreement with those 
from Refs.— ancU^. The results in this paper can be com- 
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FIG. 5: Effects of phonon bandwidth A on (a) total diffusion 
coefficient, (b) band-like contribution and (c) hopping con- 
tribution. The other control parameters are g 2 =0.5, </> 2 =0.0 
and J=0.1. 



pared with those obtained using other methods such as a 
Green's function approach^ and variational exact diago- 
nalization with a better construction of phonon states^. 
For example, with an improved numerical technique in 
Ref4i, the ID polaron mass and radius are studied for an 
expanded parameter regime, where comparisons with our 
approach are especially convenient. Despite the discrep- 
ancies between our results and those of Munn and Sil- 
bey in high temperature regime where hopping transport 
dominates, agreement is found in low temperature regime 
where band-like transport prevails. This can be under- 
stood as follows. As discussed in Reb^l, the transforma- 
tion coefficients obtained self-consistently are shown to 
have stronger wave vector dependencies as 4> increases. 
For weak off-diagonal coupling, however, the average val- 
ues of ^ and rf^ generated by the self-consistency pro- 
cedure are in reasonable agreement with those by Munn 
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FIG. 6: Comparisons of the logarithm diffusion coefficient 
of this work and those from Ref.—, where both diagonal and 
off-diagonal couplings are limited to a moderate value to guar- 
antee the definiteness of our results. 



off-diagonal coupling are shown to have qualitative agree- 
ments with those from Ref. 13 . The differences between 
the two can be ascribed to a series approximations made 
in Ref. 13 leading to the calculation of the diffusion coeffi- 
cient, which in this work are substituted with numerically 
exact evaluations. 



IV. CONCLUSIONS 

Within the framework of the Munn-Silbey approach 
to polaronic transport in organic molecular crystals, dif- 
fusion coefficients have been computed with numerical 
means while taking into account simultaneous diago- 
nal and off-diagonal exciton-phonon coupling for a wide 
range of temperatures. With the transformation coeffi- 
cients determined self-consistcntly without the additional 
approximations found in Refsj 12 ' 13 , calculated trans- 
port properties exhibit substantial corrections on those 
obtained previously thanks to the added momentum- 
dependence of the A-matrix parameters and numerically 
exact evaluation of many essential quantities. The gen- 
eral transport properties obtained are also in good qual- 
itative agreement with those obtained with the Holstein- 
Peierls model and ab-initio calculations 29 demonstrating 
robustness of the current approach. 

Comparisons between this work and the one in RefJ^ 
reveal the importance of the electron-phonon correlations 
and the polaron band structure to transport properties of 
the extended Holstein molecular crystal model. By treat- 
ing the diagonal and off-diagonal exciton-phonon cou- 
pling on an equal footing in our self-consistency proce- 
dure, the current approach is one of the few capable to 
model realistically the transport phenomenon in organic 
molecular crystals, for which the effect of off-diagonal 
coupling often can not be neglected. Off-diagonal cou- 
pling, in itself, is a transport mechanism as it is an agent 
for polaronic localization. Inclusion of off-diagonal cou- 
pling in the current approach allows us to see the critical 
role it plays in determining both polaron structures and 
transport properties over a wide range of temperatures. 
Of course, just as any other treatment of transport, the 
current approach can be improved, especially in the low 
temperature regime, where the temperature-independent 
mobility has recently been found. Work in this direction 
is in progress. 



and Silbeyi 2 - despite some systematic deviations, which 
leads to the aforementioned low-tcmpcrature agreement. 
As shown in Figs. [2] and El a higher temperature induces 
a higher amplitude of wave vector variations of the A- 
matrix parameters together with a much natter polaron 
band as compared with its Munn-Silbey counterpart, re- 
sulting in a diminished diffusion coefficient. Finally, the 
diffusion coefficients obtained here in the absence of the 
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Appendix A: Explicit expressions for the correlation functions 

The expression of the correlation functions (Vk'+q t k+qVk,k'{t)} is mainly composed of three parts as: 

/r k 1 ,k 2 V qi .q 2 (t)) = X klyk2;qit q 2 + ^Al M ;9l ,92 + A *i ,fc 2 ; 01 ,g 2 ■ (Al) 

First of all, the expression of -Xfci,fc 2 -01,92 is given by: 



-^fci,fca 591,92 



,qi;k 4 ,q 2 ) 

k 3 ,k 4 

-^y^Jk 3 UJ Q 2 f Q k 2 A ■ [A_®* (Tk 3 ,k 1 ;k 3 ,k 2 T k 4 +Q 2 ,qr,k i ,q 2 -Q 2 ) + A®^* (T k3 ,k 1 ;k 3 ,k 2 Tk 4 +Q 2 ,q 1 +Q 2 ;k 4 ,q 2 )) 
Q2 

^y!4 M Ql/ g fa ' (^-ki (Tk 3 +Q u k 1 ;k 3 M~Q 1 Tk i ,q 1 ;k 4 ,q 2 ) + A_ k * (Tfcg+Qj ,fc 1 +Q 1 ;fcs,fe 2 ^fc4,«i;fc4,?3 )) 



+ ^2 E ^i-wVl 

~i,Q 2 



x (^-fc 1 1 ^-? 1 M r fc3+Qi,fci;fe3,fe2-Qi T fc4+Q2,oi;fe4,92-Q2) + ^-fc^-g* (r , fc 3 +Q 1 ,fei;fc3,fc2-Qi T 'fc4+Q2,<ji+Q2;fe4,g2) 

+ ^ Q fc* J 4- 9 ? i 2 ( T 'fc3+Qi,fci+Qi;fc3,fc2 T fc4+Q2,9i:fc4,g2-Q2) + ^-fc*^- 91 *( T 'fc3+Qi,fci+Qi;fe3,fe2 71 fe4+Q2,gi+Q2:fc4,g2))] 5 ( A2 ) 

where 

(Tfe!, k 2 ;k 3 ,k 4 T qu ?2:?3, 94 (*)) = (^fe 1 ,fc 2 ^fe3,fe4^i,92(*)^93,94(*)) - ( 6 'fc 1 .fe 2 ^fc3,fc4)(^9 1 ,9 2 6, 93,94>: ( A3 ) 

and the expression for (0 klik2 9 k3iki 6l uq2 (t)6 q3i q 4 (t)) and (0l lM 9k 3 ,k 4 ) (#|i ,92^93, 94) nas been S iven b y Eq. [38] and 
Eq. [33l respectively. Here and onwards, AT represents the total number of sites used in the transport calculation. 

Now we turn to the expression of Y kl}k2 . >qiiq2 , which is written as: 

*Ki,fea;gi>g2 

= E {~^iy^/k 3 ^Qf Q kl (Zk 1 M;k 3 M^Q( t ))( T k 3 M\k 3 ,k 2 T k4+Q,qi\k i ,q 2 (t)) - W Ql Wq 2 

fc 3 ,fc4 ^ Q Qi,Q 2 
x (^-/n ( Z kiM+Qi\k 3 M-Qi^Q2( t ))( T k 3 +QiM\k 3 ,k 2 -Qi T k4+Q 2 ;qi;k 4 ,q 2 (t)) 

+ ^ Q fc*('^fcl+Ql,fc3+Qi;fc3,fc2V'Q2(*)>( T fc3+Ql,fcl+Qi:fc3,fe2 T fc4 + Q2,9i:fc4,92( < )>)] 

+ ~7v7 JkjUQxf-j, ■ (V'Qi^9i,fc4;fc4,92( i ))(' 7 fc3+Qi,fci;fc3,fc2 T 'fc4,9i;fc4,92(*)) - E ^Ql^Qzf-kf-L 
V Qi Qi,Q 2 

X ( j4 - g 3 1 2 (V'Ql2'9l,'i;4+Q2;fc4,92-Q2(i)>( T /C3+Ql,/c 1 ;fc3,fc2 T fe4 + Q2,9i;fe4,92-Q2( i )> 

+ ^?9* (V'Qi^9i+Q2,A;4+Q2;fc4,92(0)( r fc3+Qi,fci;fe3,fc2 r fc4+Q2,9i+Q2;fc4,92( < )))]} I ( A4 ) 

where 

(C fcl , fe2;fe3 , fe4 V9(^)> - ^{A^ fc2 [ne-^ - (1 + n)e^]S. kl+k2 , q A^^ne^ (n + l)e iw *]<5_ fc3+fe4)0 }, (A5) 



and 



(^^.faifcs^W) = ^{^^[(n + l^-ne-^.-^+fa -^: 3 fe - fe4 [(n + l)e-*-» e — ^.^^^J, (A6) 

here n is the Bose-Einstein distribution, with 2n + 1 = coth(i/3fio;). 
Finally we give the expression of ^-ki,kam,9» '■ 
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^k 1 ,k 2 \qi,q2 
k 3 ,ki Q1Q2 

x {^(-AX^^AX^'Sk^-Q^S-Q^-k.-Q, + AX%^A k _^S kl - k3 - Qu - Q2 6 Qug2 - ki 
+ Atr 2 k2 A%%^5 k2 - k3 , Q2 S- Quqi - k4 - Q2 Atr 2 k2 A k J-^S k2 . k3 , Q2 5 Quq2 . ki ) ■ [ne"*"* - (n + l)e^] 2 
+ (V'QiV'-C^W^Qi -Q2} ' ( T k 3 +QiM;k 3 M T k4+Q2,qi;k i ,q 2 )} (A7) 
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